/* Copyright (c) Colorado School of Mines, 2010.*/
/* All rights reserved.                       */

/* SUKDMIG3D: $Revision: 1.5 $ ; $Date: 2006/11/07 22:58:42 $	*/

#include "su.h"
#include "segy.h"

/*********************** self documentation **********************/
char *sdoc[]={
" 									",
"SUKDMIG3D - Kirchhoff Depth Migration of 3D poststack/prestack data	",
" 									",
" 	sukdmig3d datain= dataout= [parameters] 			",
"									",
" Required parameters:							",
" ttfile	 file for input tttables			       	",
"									",
" Optional Parameters:							",
"									",
" datain=stdin	 file for input seismic traces				",
" dataout=stdout file for common offset migration output		",
" crfile=NULL    file for cos theta and ray paths                       ",
"									",
"   The following 17 parameters describe tttables: (from ttfile header)	",
" fxgd= or from header (f1)		first x-sample in tttable	",
" nxt= or from header (ns)		number of x-samples in tttable	",
" dxgd= or from header (d1)		x-interval in tttable 		",
" fygd= or from header (f2)		first y-sample in tttable 	",
" nyt= or from header (ntr)		number of y-samples in tttable	",
" dygd= or from header (d2)		y-interval in tttable		",
" ixsf= or from header (sdel) 	        x in dxgd of first source	",
" nxs= or from header (nhs) 		number of sources in x		",
" ixsr= or from header (swdep) 	        ratio of source & gd spacing    ",
" iysf= or from header (gdel)	        y in dygd of first source 	",
" nys= or from header (nvs)		number of sources in y          ",
" iysr= or from header (gwdep)	        ratio of source & gd spacing    ",
" fzs= or from header (sdepth/1000)	first depth sample in tttable	",
" nzs= or from header (duse) 		number of depth samples in tttable",
" dzs= or from header (ep/1000)		depth interval in tttable	",
" nxgd= or from header (selev)		x size of the traveltime region ",
" nygd= or from header (gelev)		y size of the traveltime region ",
" multit= or from header (scalel)       number of multivalued traveltime",
"									",
"   The following two parameters are from data header			",
" dt= or from header (dt)	time sampling interval of input data	",
" ft= or from header (ft)	first time sample of input data		",
" dxm= or from header (d2)      mid point spacing of input data         ",
"									",
" Default: output is 5 times finer in depth and 2 times finer in lateral",
" fzo=fzs	z-coordinate of first point in output trace		",
" dzo=0.2*dzs	vertical spacing of output trace (5 times finer)	",
" nzo=5*(nzs-1)+1 number of points in output trace (5 times finer)	",	
" fxo=fxgd	x-coordinate of first output trace	 		",
" dxo=0.5*dxgd	horizontal spacing of output trace (2 times finer)	",
" nxo=2*(nxgd-1)+1 number of output traces (2 times finer)		",
" fyo=fygd	y-coordinate of first output trace			",
" dyo=0.5*dygd	horizontal spacing of output trace (2 times finer)	",
" nyo=2*(nygd-1)+1 number of output traces (2 times finer)		",
"									",
" Default: poststack migration						",	
" fxoffset=0		first offest in output in x			",
" fyoffset=0		first offest in output in y			",
" dxoffset=99999	offset increment in output in x	 		",
" dyoffset=99999	offset increment in output in y			",
" nxoffset=1		number of offsets in output in x		",
" nyoffset=1		number of offsets in output in y		",	
" xoffsetmax=99999	x-maximum absolute offset allowed in migration 	",
" yoffsetmax=99999	y-maximum absolute offset allowed in migration  ",
" xaper=nxt*dxgd/2.5	migration lateral aperature in x		",
" yaper=nyt*dygd/2.5	migration lateral aperature in y		",
" angmax=60             max angle to handle                             ",
" fmax=0.25/dt          max frequency in the data                       ",
" jpfile=stderr		job print file name 				",
" pptr=100		print verbal information at every pptr traces	",
" ntrmax=100000		maximum number of input traces to be migrated	",
" ls=0                  point =0 line source =1                         ",
"									",
" Notes:								",
" 1. Traveltime tables were generated by program SUTETRARAY (or other	",
"    ones) on very sparse tetrahedral model, with dimension nys*nxs*nzs ",
"    *nyt*nxt.                                                          ",
" 2. Input seismic traces must be SU format and can be any type of 	",
"    gathers (common shot, common offset, common CDP, and so on).	", 
" 3. Migrated traces are output in CDP gathers if velocity analysis	",
"    is required, with dimension nyoffset*nxoffset*nyo*nxo*nzo.		", 
" 4. If the offset value of an input trace is not in the offset array 	",
"    of output, the nearest one in the array is chosen. 		",
" 5. Memory requirement for this program is about			",
"    [nys*nxs*nzs*nyt*nxt+nyoffset*nxoffset*nxo*nyo*nzo+	       	",
"    nys*nxo*nzo*nyt*nxt]                                               ",
" 6. Input traces must specify source and receiver positions via header	",
"    fields tr.sx and tr.gx, as well as tr.sy and tr.gy. Offset is 	",
"    computed automatically.						",
"									",
" Disclaimer:								",
" This is a research code that will take considerable work to get into	",
" the form of a a production level 3D migration code. The code is	",
" offered as is, along with tetramod and sutetraray, to provide a starting",
" point for researchers who wish to write their own 3D migration codes.",
"									",
NULL};

/*
 * Author:  Zhaobo Meng, 01/10/97,  Colorado School of Mines 
 *
 * Trace header fields accessed: ns, dt, delrt, d2
 * Trace header fields modified: sx, gx
 */
/******************************* end self doc *******************************/

#define DEBUG

#define InfDistance 10
#define InfTime 7000 /*ms*/

#define InfByte1 255
#define InfByte2 65535

struct GD {
      float fxgd;
      float dxgd;
      int nxgd;	        /*number of samples in x in tttable*/
      float fygd;
      float dygd;
      int nygd;	        /*number of samples in tttable in y*/

      int nxt;	        /*number of lateral samples (in x) in tttable*/
      int nyt;	        /*number of lateral samples (in y) in tttable*/

      int nxo;	        /*number of samples in x in output*/
      float fxo;	/*first sample in x in output*/
      float dxo;	/*sample spacing in x in output*/
      int nyo;	        /*number of samples in y in output*/
      float fyo;	/*first sample in y in output*/
      float dyo;	/*sample spacing in y in output*/
      int nzo;	        /*number of samples in z in output*/
      float fzo;	/*first sample in z in output*/
      float dzo;	/*sample spacing in z in output*/

      float fxs;        /*first shot in x*/
      float dxs;        /*shot spacing in x*/
      int nxs;          /*number of sources in x*/
      float fys;        /*first shot in y*/
      float dys;        /*shot spacing in y*/
      int nys;          /*number of sources in y*/
      float fzs;        /*first shot in z*/
      float dzs;        /*shot spacing in z*/
      int nzs;          /*number of sources in z*/
      float xaper;	/*aperture in x*/
      float yaper;	/*aperture in y*/
      int multit;
      int ixsf;         /*first source in x in gd*/
      int ixsr;         /*source spacing/grid spacing*/
      int iysf;         /*first source in y in gd*/
      int iysr;         /*source spacing/grid spacing*/
      float sy;	        /*source position in y*/
      float gy;	        /*geophone position in y*/
      float sx;	        /*source position in x*/
      float gx;	        /*geophone position in y*/
      int nt;		/*number of samples of this trace*/
      float ft;	        /*first sample of this trace*/
      float dt;	        /*sampling rate in t*/
      float v0;
      float fmax;
      int ls;
      int ntrimax;
      FILE *crfp;
      FILE *jpfp;
};

struct TTAB {
      unsigned short t;  /*real t=t/InfByte2*InfTime*/
      unsigned char cs;  /*real cs=cs/InfByte1*/
      unsigned char r;   /*real r=r/InfByte1*InfDistance*/
};

void filt(float *trace, struct GD *gd, float *trf);

void get_bufs( struct GD *gd, float ******ttab, float *****ctab,
      float *****rtab, unsigned short ******buf, unsigned char *****cbuf,
      unsigned char  *****rbuf);

void mig3d( struct GD *gd, float *trace, float *trf, float ***mig,
      unsigned short ******buf, unsigned char  *****cbuf,
      unsigned char  *****rbuf);

void preproc( float *data, struct GD *gd);

/*segy trace*/
segy tr,tro,ttr;

int
main (int argc, char **argv)
{
      struct GD gd0;
      struct GD *gd=&gd0;

      /*************************************************************
	Memory needed for the computed area (for nxoff=nyoff=1):
		nxo * nyo * nzo * 4; (=100*100*100*4=8Mbytes)	

	Memory needed for the ttab:
		nys * nxo *nzo * nxt * nyt * 2; 
		(=20*100*200*40*40*2=1280Mbytes)

		(0,nygd)
			 ______________________________________ 
			|       |	    |		       |
			|       |	    |		       |
			|       |  ttable   |		       |
			|       |  for 1 s  |	     	       |
			|       |	    |		       |
			|_______|___________|__________________|
		(0,0) (ixsf+ixsr*ixs-nxt/2,0) (ixsf+ixsr*ixs+nxt/2,0) (nxgd,0)	

      **************************************************************/
      int nxoffset;     /*number of offsets in output in x*/
      int nyoffset;     /*number of offsets in output in y*/
      int ixoffset;     /*index for nxoffset*/
      int iyoffset;     /*index for nyoffset*/
      int ixo;	        /*index for nxo*/
      int iyo;	        /*index for nyo*/	
      int ixt;	        /*index for nxt*/
      int iyt;	        /*index for nyt*/
      int izs;	        /*index for nzs*/
      int ntrmax;	/*maximum number of input traces to be migrated*/
      int imigtr;	/*index for migrated traces*/
      int itotmigtr;    /*index for migrated traces*/
      int pptr;	        /*print verbal information for every pptr traces*/

      float fxoffset;   /*first offset in output in x*/
      float fyoffset;   /*first offset in output in y*/
      float dxm;        /*mid point spacing*/

      float angmax;     /*max angle to be able to handle*/
      float dxoffset;   /*offset increment in output in x*/
      float dyoffset;   /*offset increment in output in y*/
      float ext;	/*ending t in x=fxgd+(nxt-1)*dxgd*/
      float eyt;	/*ending t in y=fygd+(nyt-1)*dygd*/
      float ezs;	/*last time in z=fzs+(nzs-1)*dzs*/
      float ezo;	/*last point in output ezo=fzo+(nzo-1)*dzo*/
      float exo;	/*last point in output ex0=fxo+(nxo-1)*dxo*/
      float eyo;	/*last point in output ey0=fyo+(nyo-1)*dyo*/
      float exs;	/*last source exs=fxs+(nxs-1)*dxs*/
      float eys;	/*last source eys=fys+(nys-1)*dys*/
      float xoffsetmax; /*maximum absolute x offset allowed in migration*/
      float yoffsetmax; /*maximum absolute y offset allowed in migration*/
      float *****mig;   /*migrated data of size nyoffset*nxoffset*nyo*nxo*nzo*/
      float ******ttab; /*travetime table of size multit*nys*nxs*nzs*nyt*nxt*/
      float *****ctab;  /*travetime table of size nys*nxs*nzs*nyt*nxt*/
      float *****rtab;  /*travetime table of size nys*nxs*nzs*nyt*nxt*/

      unsigned short ******tbuf; /*traveltime buffer of size nys*nxo*nzo*nyt*nxt*/
      unsigned char  *****cbuf;  /*cos theta buffer of size nys*nxo*nzo*nyt*nxt*/
      unsigned char  *****rbuf;  /*cos theta buffer of size nys*nxo*nzo*nyt*nxt*/

      int iys,ixs;      /*indices for nys, nxs*/

      float *trf;       /*integrated traced*/
      int imul;

      char *datain="stdin";/*input data*/
      char *dataout="stdout";/*output data*/
      char *ttfile;     /*tttable to be read in*/
      char *crfile;     /*cos and r to be read in*/
      char *jpfile;     /*job print file name*/
      FILE *infp;	/*input file pointer*/
      FILE *outfp;	/*output file pointer*/
      FILE *ttfp;	/*tttable file pointer*/

      float xm,ym;      /*mid point*/

      /*hook up getpar to handle the parameters*/
      initargs(argc,argv);
      requestdoc(1);

      /*open input and output files*/
      if (!getparstring("datain",&datain)) {
	    infp=stdin;
      } else
	    if ((infp=fopen(datain,"r"))==NULL)
		  err("cannot open datain=%s",datain);

      if ( !getparstring("dataout",&dataout)) {
	    outfp=stdout;
      } else
	    outfp=fopen(dataout,"w");

      /***********************************************************
      Repositon a file on a logical unit: e.g. infp;
      efseeko(infp,offset,from), offset is an offset in bytes
      relative to the position specified by offset, offset=0 
      means beginning;=1 means current position;=2 means the end
      return 0 when successful, otherwise a system error code.
      ***********************************************************/
      efseeko(infp,(off_t) 0,SEEK_CUR);
      efseeko(outfp,(off_t) 0,SEEK_CUR);

      /*****************************************************
      Open tttable file, get parameters from the table
      *****************************************************/
      if ( !getparstring("ttfile",&ttfile))
	    err("must specify ttfile!");
      if ((ttfp=fopen(ttfile,"r"))==NULL)
	    err("Can not open ttfile=%s",ttfile);
      if (!fgettr(ttfp,&ttr))
	    err("Can't get first trace from tttable");

      if (getparstring("crfile",&crfile)) {
            if ((gd->crfp=fopen(crfile,"r"))==NULL)
	          err("Can not open crfile=%s",crfile);
      }
	else
	gd->crfp=NULL;
	

      if (!getparfloat("fxgd",&gd->fxgd)) gd->fxgd=ttr.f1;
      if (!getparfloat("dxgd",&gd->dxgd)) gd->dxgd=ttr.d1;
      if (!getparint("nxt",&gd->nxt)) 	  gd->nxt=ttr.ns;

      if (gd->dxgd<0.1) /* km used I guess */
            gd->v0=2.5;
      else
            gd->v0=2500.0;

      if (!getparint("multit",&gd->multit)) gd->multit=ttr.scalel;
      if (gd->multit==0) gd->multit=1;
      gd->nxt/=gd->multit;

      if (!getparint("ls",&gd->ls)) gd->ls=0;

      if (!getparfloat("fygd",&gd->fygd)) gd->fygd=ttr.f2;
      if (!getparfloat("dygd",&gd->dygd)) gd->dygd=ttr.d2;
      if (!getparint("nyt",&gd->nyt))	  gd->nyt=ttr.ntr;

      gd->nxgd=ttr.selev;
      gd->nygd=ttr.gelev;

      /*number of sources*/
      if (!getparint("nxs",&gd->nxs))	gd->nxs=ttr.nhs;
      if (!getparint("nys",&gd->nys))	gd->nys=ttr.nvs;
      if (!getparint("nzs",&gd->nzs))   gd->nzs=ttr.duse;

      /*first source*/
      if (!getparint("ixsf",&gd->ixsf))	gd->ixsf=ttr.sdel;
      if (!getparint("iysf",&gd->iysf))	gd->iysf=ttr.gdel;
      if (!getparfloat("fzs",&gd->fzs)) gd->fzs=ttr.sdepth/1000.0;

      /*source spacing*/
      if (!getparint("ixsr",&gd->ixsr))	gd->ixsr=ttr.swdep;
      if (!getparint("iysr",&gd->iysr))	gd->iysr=ttr.gwdep;
      if (!getparfloat("dzs",&gd->dzs)) gd->dzs=ttr.ep/1000.0;
     
      fclose(ttfp);
      if ((ttfp=fopen(ttfile,"r"))==NULL)
            err("Can not open ttfile=%s",ttfile);
      efseeko(ttfp,(off_t) 0,SEEK_CUR);

      if (!getparfloat("dyoffset",&dyoffset)) dyoffset=99999;

      if (!getparint("ntrmax",&ntrmax)) ntrmax=100000;
      if (!getparint("pptr",&pptr)) pptr=100;

      gd->fxs=gd->fxgd+gd->ixsf*gd->dxgd;
      gd->dxs=gd->ixsr*gd->dxgd;
      gd->fys=gd->fygd+gd->iysf*gd->dygd;
      gd->dys=gd->iysr*gd->dygd;

      exs=gd->fxs+(gd->nxs-1)*gd->dxs;	/*last source in x in tttable*/
      eys=gd->fys+(gd->nys-1)*gd->dys;	/*last source in y in tttable*/
      ezs=gd->fzs+(gd->nzs-1)*gd->dzs;	/*last source in z in tttable*/
      ext=gd->fxgd+(gd->nxgd-1)*gd->dxgd;	/*last sample in x in tttable*/
      eyt=gd->fygd+(gd->nygd-1)*gd->dygd;	/*last sample in y in tttable*/

      if (!getparint("nxo",&gd->nxo)) gd->nxo=(gd->nxgd-1)*2+1;
      if (!getparfloat("fxo",&gd->fxo)) gd->fxo=gd->fxgd;
      if (!getparfloat("dxo",&gd->dxo)) gd->dxo=gd->dxgd*0.5;

      if (!getparint("nyo",&gd->nyo)) gd->nyo=(gd->nygd-1)*2+1;
      if (!getparfloat("fyo",&gd->fyo)) gd->fyo=gd->fygd;
      if (!getparfloat("dyo",&gd->dyo)) gd->dyo=gd->dygd*0.5;

      if (!getparint("nzo",&gd->nzo)) gd->nzo=(gd->nzs-1)*5+1;
      if (!getparfloat("fzo",&gd->fzo)) gd->fzo=gd->fzs;
      if (!getparfloat("dzo",&gd->dzo)) gd->dzo=gd->dzs*0.2;

      exo=gd->fxo+(gd->nxo-1)*gd->dxo;	/*last sample in x in output*/
      eyo=gd->fyo+(gd->nyo-1)*gd->dyo;	/*last sample in y in output*/
      ezo=gd->fzo+(gd->nzo-1)*gd->dzo;	/*last sample in z in output*/

      /*****************************************************
      The migration area should lie within the traveltime
      region

	      (fxgd,fygd,fzs)
		    _________________________________
		    |(fxo,fyo,fzo)		     |
		    |       __________________       |
		    |      |	   	      |      |
		    |      |		      |      |
 		    |      |		      |      |
		    |      |		      |      |
		    |      |__________________|      |
		    |		   (exo,eyo,ezo)     |
		    |________________________________|
					       (ext,eyt,ezs)
      *****************************************************/

      /***************************************************
      The apertures for x and y is the influential region
      of the current trace
      ***************************************************/
      if (!getparfloat("xaper",&gd->xaper)) 
       	    gd->xaper=gd->nxt*gd->dxgd/2.5;
      if (!getparfloat("yaper",&gd->yaper)) 
	    gd->yaper=gd->nyt*gd->dygd/2.5;
      if (!getparfloat("xoffsetmax",&xoffsetmax)) xoffsetmax=3000;
      if (!getparfloat("yoffsetmax",&yoffsetmax)) yoffsetmax=3000;

      /**************************************************
      For effective MVA, we need some nxoffset, nyoffset 
      about 10
      **************************************************/
      if (!getparint("nxoffset",&nxoffset)) nxoffset=1;
      if (!getparint("nyoffset",&nyoffset)) nyoffset=1;
      if (!getparfloat("fxoffset",&fxoffset)) fxoffset=0.;
      if (!getparfloat("fyoffset",&fyoffset)) fyoffset=0.;
      if (!getparfloat("dxoffset",&dxoffset)) dxoffset=99999;

      /*********************************************************
      Get information from the first header of data
      *********************************************************/
      if (!fgettr(infp,&tr)) 
	    err("Can't get first trace");
      gd->nt=tr.ns;	/*number of samples*/
      if (!getparfloat("dt",&gd->dt)) 
	    gd->dt=((double)tr.dt)/1000.0;  /*ms*/
      if (gd->dt<0.0000001) 
	    err("dt must be positive!");
      if (!getparfloat("ft",&gd->ft)) 
	    gd->ft=tr.delrt/1000.0;
      if (!getparfloat("dxm",&dxm))
            dxm=tr.d2;
      if (dxm<0.0000001) err("dxm must be positive");
      if (!getparfloat("angmax",&angmax))
            angmax=60.0*PI/180.0;
      if (angmax<0.0000001) err("angmax must be positive");
      gd->ntrimax=2*dxm*sin(angmax)/(gd->v0*gd->dt);
      if (gd->ntrimax<1) gd->ntrimax=1;

      if (!getparfloat("fmax",&gd->fmax)) gd->fmax=0.25/gd->dt;

      /*********************************************************
      Preparing the job print
      *********************************************************/
      if ( !getparstring("jpfile",&jpfile)) {
	    gd->jpfp=stderr;
      } else
	    gd->jpfp=fopen(jpfile,"w");

      fprintf(gd->jpfp,"\n");
      fprintf(gd->jpfp," Migration parameters\n");
      fprintf(gd->jpfp," ================\n");
      fprintf(gd->jpfp," ttfile=%s \n",ttfile);
      fprintf(gd->jpfp," \n");
      fprintf(gd->jpfp," nxt=%d fxgd=%g dxgd=%g\n",gd->nxt,gd->fxgd,gd->dxgd);
      fprintf(gd->jpfp," nyt=%d fygd=%g dygd=%g\n",gd->nyt,gd->fygd,gd->dygd);
      fprintf(gd->jpfp," nxs=%d fxs=%g dxs=%g\n",gd->nxs,gd->fxs,gd->dxs);
      fprintf(gd->jpfp," nys=%d fys=%g dys=%g\n",gd->nys,gd->fys,gd->dys);
      fprintf(gd->jpfp," nzs=%d gd->fzs=%g dzs=%g\n",gd->nzs,gd->fzs,gd->dzs);
      fprintf(gd->jpfp," nxgd=%d nygd=%d\n\n",gd->nxgd,gd->nygd);
      fprintf(gd->jpfp," nzo=%d fzo=%g dzo=%g\n",gd->nzo,gd->fzo,gd->dzo);
      fprintf(gd->jpfp," nxo=%d fxo=%g dxo=%g\n",gd->nxo,gd->fxo,gd->dxo);
      fprintf(gd->jpfp," nyo=%d fyo=%g dyo=%g\n\n",gd->nyo,gd->fyo,gd->dyo);
      fprintf(gd->jpfp," nt=%d ft=%g dt=%g dxm=%g\n\n",gd->nt,gd->ft,gd->dt,dxm);
      fprintf(gd->jpfp," nxoffset=%d fxoffset=%g dxoffset=%g\n",
	    nxoffset,fxoffset,dxoffset);
      fprintf(gd->jpfp," nyoffset=%d fyoffset=%g dyoffset=%g\n",
	    nyoffset,fyoffset,dyoffset);
      fprintf(gd->jpfp," xaper=%g xoffsetmax=%g\n",
	    gd->xaper,xoffsetmax);
      fprintf(gd->jpfp," yaper=%g yoffsetmax=%g\n",
	    gd->yaper,yoffsetmax);
      fprintf(gd->jpfp," ntrmax=%d pptr=%d\n",ntrmax,pptr);
      fprintf(gd->jpfp," ext=%g exo=%g\n",ext,exo);
      fprintf(gd->jpfp," eyt=%g eyo=%g\n",eyt,eyo);
      fprintf(gd->jpfp," ezs=%g ezo=%g\n",ezs,ezo);
      fprintf(gd->jpfp," ====================================\n");
      fflush(gd->jpfp);

      if (gd->fxgd>gd->fxo || ext<exo ||
	  gd->fygd>gd->fyo || eyt<eyo) {
            fprintf(gd->jpfp,"fxgd=%g > fxo=%g ?\n",gd->fxgd,gd->fxo);
            fprintf(gd->jpfp,"ext=%g < exo=%g ?\n",ext,exo);
            fprintf(gd->jpfp,"fygd=%g > fyo=%g ?\n",gd->fygd,gd->fyo);
            fprintf(gd->jpfp,"eyt=%g < eyo=%g ?\n",eyt,eyo);
	    err(" migration output range is out of tttable");
      }

      /****************************************************
      Allocate space
      ****************************************************/
      fprintf(gd->jpfp,"allocating mig of size %g Mbytes\n",(4.0e-6)*
	    gd->nzo*gd->nxo*gd->nyo*nxoffset*nyoffset);

      mig =ealloc5float(gd->nzo,gd->nxo,gd->nyo,nxoffset,nyoffset);

      fprintf(gd->jpfp,"allocating ttab, ctab and rtab of size %g Mbytes\n",(12.0e-6)*
	    gd->nzs*gd->nxs*gd->nys*gd->nxt*gd->nyt);
      ttab=ealloc6float(gd->nxt,gd->nyt,gd->nzs,gd->nxs,gd->nys,gd->multit);
      ctab=ealloc5float(gd->nxt,gd->nyt,gd->nzs,gd->nxs,gd->nys);
      rtab=ealloc5float(gd->nxt,gd->nyt,gd->nzs,gd->nxs,gd->nys);

      fprintf(gd->jpfp,"allocating tbuf, cbuf and rbuf of size %g Mbytes\n",(4.0e-6)*
	    gd->nys*gd->nxo*gd->nzo*gd->nyt*gd->nxt);
      tbuf=ealloc6ushort(gd->multit,gd->nxt,gd->nyt,gd->nzo,gd->nxo,gd->nys);
      cbuf=ealloc5uchar(gd->nxt,gd->nyt,gd->nzo,gd->nxo,gd->nys);
      rbuf=ealloc5uchar(gd->nxt,gd->nyt,gd->nzo,gd->nxo,gd->nys);

      trf=ealloc1float(gd->nt + 2*gd->ntrimax);

      memset((void *)mig[0][0][0][0],0,
	    nxoffset*nyoffset*gd->nxo*gd->nyo*gd->nzo*sizeof(float));

      fprintf(gd->jpfp,"\nReading ttables ...");
      fprintf(gd->jpfp,"\n");

      /*****************************************************		       
      Compute traveltime residual, which is the difference
      of the tttable and the traveltimes calculated
      from the reference model	
      *****************************************************/
      for (iys=0;iys<gd->nys;iys++) {
	    for(ixs=0;ixs<gd->nxs;ixs++){			
		  for (izs=0;izs<gd->nzs;izs++) {
			fprintf(gd->jpfp,
			      "Reading source No. (%d,%d,%d)\n",ixs,iys,izs);
			for (iyt=0;iyt<gd->nyt;iyt++) {
			      if (!fgettr(ttfp,&ttr))
				    err("Can't read ttfile");
			      for (imul=0;imul<gd->multit;imul++) {
			            for (ixt=0;ixt<gd->nxt;ixt++) {
				          ttab[imul][iys][ixs][izs][iyt][ixt]
					        =ttr.data[ixt+imul*gd->nxt];
			            }
			      }
                              if (gd->crfp!=NULL) {
                                    if (fread(ctab[iys][ixs][izs][iyt],
					  sizeof(float),gd->nxt,gd->crfp)!=gd->nxt) 
				          err("Can not read in crfile");
				    for (ixt=0;ixt<gd->nxt;ixt++) 
				          ctab[iys][ixs][izs][iyt][ixt]=1.0;

				    if (fread(rtab[iys][ixs][izs][iyt],
					  sizeof(float),gd->nxt,gd->crfp)!=gd->nxt) 
				          err("Can not read in crfile");
				    for (ixt=0;ixt<gd->nxt;ixt++) 
				          rtab[iys][ixs][izs][iyt][ixt]=1.0;
			      }
			}
		  }
	    }
      }

      fclose(ttfp);
      fprintf(gd->jpfp,"\nTable reading done\n");	

      fprintf(gd->jpfp,"\nGet traveltime buffer\n");
      get_bufs(
            gd, 
            ttab,       /*tttable, of size multit*nys*nxs*nzs*nyt*gd->nxt*/
            ctab,       /*tttable, of size nys*nxs*nzs*nyt*nxt*/
            rtab,       /*tttable, of size nys*nxs*nzs*nyt*nxt*/
            tbuf,       /*shot at nys*nxo*nzo*nyt*nxt*/
	    cbuf,       /*shot at nys*nxo*nzo*nyt*nxt*/
            rbuf);      /*shot at nys*nxo*nzo*nyt*nxt*/

      fprintf(gd->jpfp,"Start migration ...\n");
      fflush(gd->jpfp);
	
      imigtr=1;
      itotmigtr=0; 

      do { 
	    /*****************************************************	
	    To determine offset index
	    *****************************************************/
	    gd->sx=tr.sx;    /*source position in x*/
	    gd->gx=tr.gx;    /*geophone position in x*/
	    gd->sy=tr.sy;    /*source position in y*/
	    gd->gy=tr.gy;    /*geophone position in y*/ 

            #ifdef DEBUG
	    gd->sy=3;	
	    gd->gy=3;
	    fprintf(gd->jpfp,"trace %d: sx,gx,sy,gy=%g %g %g %g\n",
		  imigtr,gd->sx,gd->gx,gd->sy,gd->gy);
	    fprintf(gd->jpfp,"trace in line: %d; trace in reel: %d\n",
		  tr.tracl,tr.tracr);
            #endif

            xm=0.5*(gd->sx+gd->gx);
            ym=0.5*(gd->sy+gd->gy);

	    if (xm<gd->fxs || xm>exs || 
	        MAX(gd->gx-gd->sx,gd->sx-gd->gx)>xoffsetmax  || 
		ym<gd->fys || ym>eys || 
		MAX(gd->gy-gd->sy,gd->sy-gd->gy)>yoffsetmax )
		  continue; 

	    /****************************************************
	    Determine the index for offset in x, fxoffset is the 
	    starting offset
	    ****************************************************/
	    ixoffset=(int)((gd->gx-gd->sx-fxoffset)/dxoffset+0.5);
	    if (ixoffset<0) ixoffset=0;	
	    if (ixoffset>=nxoffset) ixoffset=nxoffset-1;

	    /****************************************************
	    Determine the index for offset in y, fxoffset is the 
	    starting offset
	    ****************************************************/
	    iyoffset=(int)((gd->gy-gd->sy-fyoffset)/dyoffset+0.5);
	    if (iyoffset<0) iyoffset=0;	
	    if (iyoffset>=nyoffset) iyoffset=nyoffset-1;

	    /*******************************************
	    Take the t-derivative of the data
	    *******************************************/
	    /*For the 1D synthetic data, turn this off*/

            filt(tr.data,gd,trf); 

	    mig3d(gd,
	       	  tr.data,  /*the data of this trace*/
                  trf,      /*integrated trace*/
		  mig[iyoffset][ixoffset],/*of nyo*nxo*nzo*/
		  tbuf,	    
		  cbuf,
		  rbuf); 

	    itotmigtr++;
	    if ((imigtr-1)%pptr ==0 ){
		  fprintf(gd->jpfp,"Migrated trace %d\n",imigtr);
		  fflush(gd->jpfp);
    	    }
	    imigtr++;
      } while (fgettr(infp,&tr) && imigtr<ntrmax);

      fprintf(gd->jpfp,"Migrated %d traces in total\n",itotmigtr);

      memset((void *)&tro,0,sizeof(segy));
      tro.ns=gd->nzo;
      tro.dt=1000*(int)(1000*gd->dt+0.5);
      tro.delrt=0.0;
      tro.f1=gd->fzo;
      tro.d1=gd->dzo;
      tro.f2=gd->fxo;
      tro.d2=gd->dxo;
      for (iyo=0;iyo<gd->nyo;iyo++) {
	    for (ixo=0;ixo<gd->nxo;ixo++) {
		  for(iyoffset=0;iyoffset<nyoffset;iyoffset++) {
			for(ixoffset=0;ixoffset<nxoffset;ixoffset++)  {
			      memcpy((void *)tro.data,(const void*)
			      mig[iyoffset][ixoffset][iyo][ixo],
			      gd->nzo*sizeof(float));
			      tro.offset=fxoffset+ixoffset*dxoffset;
			      tro.tracr=1+ixo;
			      tro.tracl=1+ixoffset+nxoffset*ixo;
			      tro.cdp=gd->fxo+ixo*gd->dxo;
			      tro.cdpt=1+ixoffset;

			      /*write out*/
			      fputtr(outfp,&tro);
			}
		  }
	    }
      }

      fprintf(gd->jpfp,"\n");
      fprintf(gd->jpfp,"Output done\n");
      fflush(gd->jpfp);

      efclose(gd->jpfp);
      efclose(outfp);

      free6ushort(tbuf);
      free5uchar(cbuf);
      free5uchar(rbuf);
      free6float(ttab);
      free5float(ctab);
      free5float(rtab);
      free5float(mig);
      return(CWP_Exit());
}

void mig3d( struct GD *gd, float *trace, float *trf, float ***mig,
      unsigned short ******tbuf, unsigned char  *****cbuf,
      unsigned char  *****rbuf)
/******************************************************************
mig3d - migrate a single trace
*******************************************************************
Function prototype:
void mig3d( struct GD *gd, float *trace, float *trf, float ***mig,
      unsigned short ******tbuf, unsigned char  *****cbuf,
      unsigned char  *****rbuf);
*******************************************************************
Inputs:
struct *gd   grid information
float *trace trace data
float *trf   integrated trace
unsigned short *******tbuf traveltime buffer
unsigned char  *****cbuf  cosine buffer
unsigned char  *****rbuf  path buffer
Outputs:
float ***mig migrated section

*******************************************************************
Notes:

For formulations:
		    |    ps  + pg    |
	h(y,xi)=det | d(ps+pg)/d xi1 |
		    | d(ps+pg)/d xi2 |

	       = 2 cos(theta)^2 z (rs+rg)(rs^2+rg^2)/(c^3 rs^3 rg^3);

	a(y,xi)=A(x,xs)A(x,xg)
	       =1/(16 pi^2 rs rg);

	|\Delta_y phi|= 2cos(theta)/c(y);

	here y is (y_1,y_2,y_3) meaning the outputing points.

	The correct correction factor is:

	| h |/a /|\Delta_y phi| =
	16 pi^2 cos(theta) y_3 (r_s+rg)(rs^2+rg^2)/(v^2 rs^2 rg^2); 

*******************************************************************
Author: CWP: Zhaobo Meng, Sept 1997
*******************************************************************/

{
      int nyof;	        /*starting index for y scope: 0<=nyof<nyo*/
      int nyoe;	        /*ending index for y scope: 0<=nyoe<nyo*/
      int nxof;	        /*starting index for x scope: 0<=nxof<nx */
      int nxoe;	        /*ending index for x scope: 0<=nxoe<nx */
      int iyo;	        /*index for nyof and nyoe*/
      int ixo;	        /*index for nxof and nxoe*/
      float ym;	        /*y midpoint position*/
      float xm;	        /*x midpoint position*/
      float yo; 	/*y position */
      float xo;	        /*x position */

      float xs,ys,ybeta,ybeta0,xalpha,xalpha0;
      int ixs,iys,izo;
      float xt,yt;
      float xs0a,xs0a0,xg0a,xg0a0,xs1a,xs1a0,xg1a,xg1a0;
      float ys0b,ys0b0,yg0b,yg0b0,ys1b,ys1b0,yg1b,yg1b0;
      float tiys0,tiys1,ttotal=0.0;
      float tiyg0,tiyg1;
      float ciys0,ciys1;
      float riys0,riys1,riyg0,riyg1,rtotal;
      float cosine;
      int itintr;
      float tintr,wtintr,wtintr0;
      float amp=1.0;
      int ixts0,ixtg0,ixts1,ixtg1;
      int iyts0,iytg0,iyts1,iytg1;
      float rs,rg;
      float zo;
      int im;

      /*for the meaning of these varibles, refer to Claerbout's
      triabgle filter notation*/
      float **rhos,**rhog;
      int ntri;    /*triangle filter length*/
      float drho,ddxg,ddyg,ddxs,ddys;
      float temp,ts=0.0,tg=0.0;

      rhos=ealloc2float(gd->nxo,gd->nyo);
      rhog=ealloc2float(gd->nxo,gd->nyo);

      /**************************************************************
      Geometry:

		fz _______________________
		   |			 |
		   |fzs_______________   |
		   |   |	     |   |
		   |   |	     |   |
		   |   |   Ttimes    |   |
		   |   |	     |   |	
		   |   |_____________|   |
		   |			 |
		   |			 |
		   -----------------------
	
      and the y dimension is similar to the x dimension
      **************************************************************/

      ym=0.5*(gd->sy+gd->gy);		/*mid point in y*/
      xm=0.5*(gd->sx+gd->gx);		/*mid point in x*/

      #ifdef DEBUG
      fprintf(gd->jpfp,"xm,ym=%g %g\n",xm,ym);
      #endif

      /*******************************************************
      For all ixo from nxof to nxoe and iyo from nyof to nyoe:
      ON OUTPUT GRIDS
      ********************************************************/
      nyof=(ym-gd->yaper-gd->fyo)/gd->dyo+0.5;
      if (nyof<0) nyof=0;
      nyoe=(ym+gd->yaper-gd->fyo)/gd->dyo+0.5;
      if (nyoe>=gd->nyo) nyoe=gd->nyo-1;

      #ifdef DEBUG
      fprintf(gd->jpfp,"This trace: fyo,dyo,nyof,nyoe=%g %g %d %d\n",
      gd->fyo,gd->dyo,nyof,nyoe);
      #endif

      nxof=(xm-gd->xaper-gd->fxo)/gd->dxo+0.5;
      if (nxof<0) nxof=0;
      nxoe=(xm+gd->xaper-gd->fxo)/gd->dxo+0.5;
      if (nxoe>=gd->nxo) nxoe=gd->nxo-1;

      #ifdef DEBUG
      fprintf(gd->jpfp,"This trace: fxo,dxo,nxof,nxoe=%g %g %d %d\n",
      gd->fxo,gd->dxo,nxof,nxoe);
      #endif

      for (iyo=nyof;iyo<=nyoe;iyo++){	
	    yo=gd->fyo+iyo*gd->dyo;
    	    ys=(yo-gd->fys)/gd->dys;
    	    iys=(int)ys;
	    iys=MIN(iys,gd->nys-2);
	    iys=MAX(0,iys);
	    ybeta=ys-iys;
            ybeta0=1.0-ybeta;

            yt=(gd->sy-gd->fygd-ybeta*gd->dys)/gd->dygd;
            iyts0=(int)yt;
            iyts0=MIN(iyts0,gd->nygd-2);
            iyts0=MAX(iyts0,0);
            ys0b=yt-iyts0;
            ys0b0=1.0-ys0b;
	    iyts0-=(gd->iysf+gd->iysr*iys-gd->nyt/2);
	    iyts0=MIN(iyts0,gd->nyt-2);
            iyts0=MAX(iyts0,0);

            yt=(gd->gy-gd->fygd-ybeta*gd->dys)/gd->dygd;
            iytg0=(int)yt;
            iytg0=MIN(iytg0,gd->nygd-2);
            iytg0=MAX(iytg0,0);
            yg0b=yt-iytg0;
            yg0b0=1.0-yg0b;
	    iytg0-=(gd->iysf+gd->iysr*iys-gd->nyt/2);
	    iytg0=MIN(iytg0,gd->nyt-2);
            iytg0=MAX(iytg0,0);

            yt=(gd->sy-gd->fygd+ybeta0*gd->dys)/gd->dygd;
            iyts1=(int)yt;
            iyts1=MIN(iyts1,gd->nygd-2);
            iyts1=MAX(iyts1,0);
            ys1b=yt-iyts1;
            ys1b0=1.0-ys1b;
	    iyts1-=(gd->iysf+gd->iysr*(iys+1)-gd->nyt/2);
	    iyts1=MIN(iyts1,gd->nyt-2);
            iyts1=MAX(iyts1,0);

            yt=(gd->gy-gd->fygd+ybeta0*gd->dys)/gd->dygd;
            iytg1=(int)yt;
            iytg1=MIN(iytg1,gd->nygd-2);
            iytg1=MAX(iytg1,0);
            yg1b=yt-iytg1;
            yg1b0=1.0-yg1b;
	    iytg1-=(gd->iysf+gd->iysr*(iys+1)-gd->nyt/2);
	    iytg1=MIN(iytg1,gd->nyt-2);
            iytg1=MAX(iytg1,0);

	    for (ixo=nxof;ixo<=nxoe;ixo++){
		  xo=gd->fxo+ixo*gd->dxo;
    	          xs=(xo-gd->fxs)/gd->dxs;
    	          ixs=(int)xs;
	          ixs=MIN(ixs,gd->nxs-2);
	          ixs=MAX(0,ixs);
		  xalpha=xs-ixs;
                  xalpha0=1.0-xalpha;

                  /*to understand this following paragraph, see
                   Claerbout's notation*/

                  rhos[iyo][ixo]=sqrt((xo-gd->sx)*(xo-gd->sx)+(yo-gd->sy)*(yo-gd->sy));
                  rhog[iyo][ixo]=sqrt((xo-gd->gx)*(xo-gd->gx)+(yo-gd->gy)*(yo-gd->gy));

                  ddxs=fabs(xo-gd->sx)*gd->dxo/rhos[iyo][ixo];
                  ddys=fabs(yo-gd->sy)*gd->dyo/rhog[iyo][ixo];

                  ddxg=fabs(xo-gd->gx)*gd->dxo/rhos[iyo][ixo];
                  ddyg=fabs(yo-gd->gy)*gd->dyo/rhog[iyo][ixo];

                  drho=sqrt(0.5*(ddxs*ddxs+ddys*ddys+ddxg*ddxg+ddyg*ddyg));
                  temp=4.0/gd->v0/gd->v0*drho/gd->dt;
                  rhos[iyo][ixo]*=temp;
                  rhog[iyo][ixo]*=temp;

                  xt=(gd->sx-gd->fxgd-xalpha*gd->dxs)/gd->dxgd;
                  ixts0=(int)xt;
	      	  ixts0=MIN(ixts0,gd->nxgd-2);
                  ixts0=MAX(ixts0,0);
                  xs0a=xt-ixts0;
                  xs0a0=1.0-xs0a;
                  ixts0-=(gd->ixsf+gd->ixsr*ixs-gd->nxt/2);
		  ixts0=MIN(ixts0,gd->nxt-2);
                  ixts0=MAX(ixts0,0);

                  xt=(gd->gx-gd->fxgd-xalpha*gd->dxs)/gd->dxgd;
                  ixtg0=(int)xt;
		  ixtg0=MIN(ixtg0,gd->nxgd-2);
                  ixtg0=MAX(ixtg0,0);
                  xg0a=xt-ixtg0;
                  xg0a0=1.0-xg0a;
                  ixtg0-=(gd->ixsf+gd->ixsr*ixs-gd->nxt/2);
		  ixtg0=MIN(ixtg0,gd->nxt-2);
                  ixtg0=MAX(ixtg0,0);

                  xt=(gd->sx-gd->fxgd+xalpha0*gd->dxs)/gd->dxgd;
                  ixts1=(int)xt;
		  ixts1=MIN(ixts1,gd->nxgd-2);
                  ixts1=MAX(ixts1,0);
                  xs1a=xt-ixts1;
                  xs1a0=1.0-xs1a;
                  ixts1-=(gd->ixsf+gd->ixsr*(ixs+1)-gd->nxt/2);
		  ixts1=MIN(ixts1,gd->nxt-2);
                  ixts1=MAX(ixts1,0);

                  xt=(gd->gx-gd->fxgd+xalpha0*gd->dxs)/gd->dxgd;
                  ixtg1=(int)xt;
		  ixtg1=MIN(ixtg1,gd->nxgd-2);
                  ixtg1=MAX(ixtg1,0);
                  xg1a=xt-ixtg1;
                  xg1a0=1.0-xg1a;
                  ixtg1-=(gd->ixsf+gd->ixsr*(ixs+1)-gd->nxt/2);
		  ixtg1=MIN(ixtg1,gd->nxt-2);
                  ixtg1=MAX(ixtg1,0);

		  /********************************
		  Apply anti-aliasing operator here
		  using scaling functions
		  ********************************/
                  for(izo=0;izo<gd->nzo;izo++) {
                        ntri=MAX(1,(rhos[iyo][ixo]/ts+rhog[iyo][ixo]/tg)-1);

                        zo=gd->fzo+izo*gd->dzo;

                        if (gd->crfp!=NULL) {
                              /*cosine table at [iys][ixo][izo][nyt][nxt]*/
                              ciys0=yg0b0*(xg0a0*cbuf[iys][ixo][izo][iytg0][ixtg0]+
                                         xg0a *cbuf[iys][ixo][izo][iytg0][ixtg0+1])+
			            yg0b *(xg0a0*cbuf[iys][ixo][izo][iytg0+1][ixtg0]+
                                         xg0a *cbuf[iys][ixo][izo][iytg0+1][ixtg0+1]);

                              /*cos table at [iys+1][ixo][izo][nyt][nxt]*/
                              ciys1=yg1b0*(xg1a0*cbuf[iys+1][ixo][izo][iytg1][ixtg1]+
                                          xg1a *cbuf[iys+1][ixo][izo][iytg1][ixtg1+1])+
			             yg1b *(xg1a0*cbuf[iys+1][ixo][izo][iytg1+1][ixtg1]+
                                          xg1a *cbuf[iys+1][ixo][izo][iytg1+1][ixtg1+1]);

			      cosine=(ybeta0*ciys0+ybeta*ciys1)/(float)InfByte1;

                              /*r table at [iys][ixo][izo][nyt][nxt]*/
                              riys0=ys0b0*(xs0a0*rbuf[iys][ixo][izo][iyts0][ixts0]+
                                         xs0a *rbuf[iys][ixo][izo][iyts0][ixts0+1])+
			            ys0b *(xs0a0*rbuf[iys][ixo][izo][iyts0+1][ixts0]+
                                         xs0a *rbuf[iys][ixo][izo][iyts0+1][ixts0+1]);

                              riyg0=yg0b0*(xg0a0*rbuf[iys][ixo][izo][iytg0][ixtg0]+
                                         xg0a *rbuf[iys][ixo][izo][iytg0][ixtg0+1])+
			            yg0b *(xg0a0*rbuf[iys][ixo][izo][iytg0+1][ixtg0]+
                                         xg0a *rbuf[iys][ixo][izo][iytg0+1][ixtg0+1]);

                              /*traveltime table at [iys+1][ixo][izo][nyt][nxt]*/
                              riys1=ys1b0*(xs1a0*rbuf[iys+1][ixo][izo][iyts1][ixts1]+
                                         xs1a *rbuf[iys+1][ixo][izo][iyts1][ixts1+1])+
			            ys1b *(xs1a0*rbuf[iys+1][ixo][izo][iyts1+1][ixts1]+
                                         xs1a *rbuf[iys+1][ixo][izo][iyts1+1][ixts1+1]);

                              riyg1=yg1b0*(xg1a0*rbuf[iys+1][ixo][izo][iytg1][ixtg1]+
                                          xg1a *rbuf[iys+1][ixo][izo][iytg1][ixtg1+1])+
			             yg1b *(xg1a0*rbuf[iys+1][ixo][izo][iytg1+1][ixtg1]+
                                          xg1a *rbuf[iys+1][ixo][izo][iytg1+1][ixtg1+1]);
                              
			      rs=(ybeta0*riys0+ybeta*riys1)/
				     (float)InfByte1*InfDistance;

			      rg=(ybeta0*riyg0+ybeta*riyg1)/
				     (float)InfByte1*InfDistance;

                              rtotal=rs+rg;
                        } else {
			      cosine=1.0;

                              rs=sqrt((xo-gd->sx)*(xo-gd->sx)+
                                      (yo-gd->sy)*(yo-gd->sy)+
                                      zo*zo);

                              rg=sqrt((xo-gd->gx)*(xo-gd->gx)+
                                      (yo-gd->gy)*(yo-gd->gy)+
                                      zo*zo);

			      rtotal=rs+rg;

                              #ifdef DEBUG
                              if (fabs((rs+rg)*1000.0-ttotal)/((rs+rg)*1000.0)>0.05 
				    && izo==0)
                              fprintf(gd->jpfp,
				    "at: ixo=%d,iyo=%d,izo=%d: est. ttot=%f, err=%f\n",
                              ixo,iyo,izo,
                              ttotal,((rs+rg)*1000.0-ttotal)/((rs+rg)*1000.0));
                              #endif
                        }

			/*cos(theta) y_3 (r_s+rg)(rs^2+rg^2)/(v^2 rs^2 rg^2)*/
                        amp=cosine*zo*rtotal*(1.0/rs/rs+1.0/rg/rg);

                        for (im=0;im<gd->multit;im++) {
                              if (tbuf[im][iys][ixo][izo][iyts0][ixts0]>=InfByte2     ||
		                  tbuf[im][iys][ixo][izo][iyts0+1][ixts0]>=InfByte2   ||
		                  tbuf[im][iys][ixo][izo][iyts0][ixts0+1]>=InfByte2   ||
		                  tbuf[im][iys][ixo][izo][iyts0+1][ixts0+1]>=InfByte2 ||
		                  tbuf[im][iys][ixo][izo][iytg0][ixtg0]>=InfByte2     ||
		                  tbuf[im][iys][ixo][izo][iytg0+1][ixtg0]>=InfByte2   ||
		                  tbuf[im][iys][ixo][izo][iytg0][ixtg0+1]>=InfByte2   ||
		                  tbuf[im][iys][ixo][izo][iytg0+1][ixtg0+1]>=InfByte2 ||
		                  tbuf[im][iys+1][ixo][izo][iyts1][ixts1]>=InfByte2   ||
		                  tbuf[im][iys+1][ixo][izo][iyts1+1][ixts1]>=InfByte2 ||
		                  tbuf[im][iys+1][ixo][izo][iyts1][ixts1+1]>=InfByte2 ||
		                  tbuf[im][iys+1][ixo][izo][iyts1+1][ixts1+1]>=InfByte2 ||
		                  tbuf[im][iys+1][ixo][izo][iytg1][ixtg1]>=InfByte2   ||
		                  tbuf[im][iys+1][ixo][izo][iytg1+1][ixtg1]>=InfByte2 ||
		                  tbuf[im][iys+1][ixo][izo][iytg1][ixtg1+1]>=InfByte2 ||
		                  tbuf[im][iys+1][ixo][izo][iytg1+1][ixtg1+1]>=InfByte2) {
                                    fprintf(gd->jpfp," skipped ixo=%d ",ixo);
                                    continue;
                              }

                              /*traveltime table at [iys][ixo][izo][nyt][nxt]*/
                              tiys0=ys0b0*(xs0a0*tbuf[im][iys][ixo][izo][iyts0][ixts0]+
                                           xs0a *tbuf[im][iys][ixo][izo][iyts0][ixts0+1])+
			            ys0b *(xs0a0*tbuf[im][iys][ixo][izo][iyts0+1][ixts0]+
                                           xs0a *tbuf[im][iys][ixo][izo][iyts0+1][ixts0+1]);

                              tiyg0=yg0b0*(xg0a0*tbuf[im][iys][ixo][izo][iytg0][ixtg0]+
                                           xg0a *tbuf[im][iys][ixo][izo][iytg0][ixtg0+1])+
			            yg0b *(xg0a0*tbuf[im][iys][ixo][izo][iytg0+1][ixtg0]+
                                           xg0a *tbuf[im][iys][ixo][izo][iytg0+1][ixtg0+1]);

                              /*traveltime table at [iys+1][ixo][izo][nyt][nxt]*/
                              tiys1=ys1b0*(xs1a0*tbuf[im][iys+1][ixo][izo][iyts1][ixts1]+
                                           xs1a *tbuf[im][iys+1][ixo][izo][iyts1][ixts1+1])+
			            ys1b *(xs1a0*tbuf[im][iys+1][ixo][izo][iyts1+1][ixts1]+
                                       xs1a *tbuf[im][iys+1][ixo][izo][iyts1+1][ixts1+1]);

                              tiyg1=yg1b0*(xg1a0*tbuf[im][iys+1][ixo][izo][iytg1][ixtg1]+
                                           xg1a *tbuf[im][iys+1][ixo][izo][iytg1][ixtg1+1])+
			            yg1b *(xg1a0*tbuf[im][iys+1][ixo][izo][iytg1+1][ixtg1]+
                                        xg1a *tbuf[im][iys+1][ixo][izo][iytg1+1][ixtg1+1]);
                              ts=(ybeta0*tiys0+ybeta*tiys1)/(float)InfByte2*InfTime;
                              tg=(ybeta0*tiyg0+ybeta*tiyg1)/(float)InfByte2*InfTime;
                        
                              ts=MAX(ts,0.001);
                              tg=MAX(tg,0.001);

                              ttotal=ts+tg;

                              if (ttotal>=InfTime ||
                                    ttotal<=0.0) continue;

	                      tintr=(ttotal-gd->ft)/gd->dt;
	                      itintr=(int)tintr;

	                      if (itintr-ntri>0 && itintr+ntri+1 < gd->nt){
	       	                    wtintr=tintr-itintr;
		                    wtintr0=1.0-wtintr;

 		                    mig[iyo][ixo][izo]+=amp*(wtintr0*
                           (-trf[itintr-ntri]+2.0*trf[itintr]-trf[itintr+ntri])
		                          +wtintr *
                           (-trf[itintr-ntri+1]+2.0*trf[itintr+1]-trf[itintr+ntri+1]));
			      }
                        }
                  }
	    }
      }
}

void get_bufs( struct GD *gd, float ******ttab, float *****ctab,
      float *****rtab, unsigned short ******tbuf, unsigned char *****cbuf,
      unsigned char  *****rbuf)
/**********************************************************************
get_bufs - Traveltime buffer. This will dramatically speed up the code.
If you could even give a larger buffer (use nyo instead of
nys), you can further speed up the code.
***********************************************************************
Function prototype:
void get_bufs( struct GD *gd, float ******ttab, float *****ctab,
      float *****rtab, unsigned short ******tbuf, unsigned char *****cbuf,
      unsigned char  *****rbuf);
***********************************************************************
Inputs:
struct *gd   grid information
float ******ttab multivalued traveltime table

outputs:
float *****ctab           cosine buffer
float *****rtab           ray path buffer
unsigned short ******tbuf multivalued traveltime buffer
***********************************************************************
Author: CWP: Zhaobo Meng, Sept 1997
***********************************************************************/
{
      int iyt;	        /*index for nyt*/
      int ixt;	        /*index for nxt*/
      int iys;          /*index for nys*/
      int ixs;          /*index for nxs*/
      int izs;	        /*index for nzs*/
      int izo;          /*index for nzo*/
      int ixo;          /*index for nxo*/      
      int im;

      float xalpha,xalpha0;
      float zgamma,zgamma0;	

      float xo,xs;
      float zo,zs;
      float t,c,r;

      for (ixo=0;ixo<gd->nxo;ixo++) {
            xo=gd->fxo+ixo*gd->dxo;
    	    xs=(xo-gd->fxs)/gd->dxs;
            ixs=(int)xs;
            ixs=MIN(ixs,gd->nxs-2);
            ixs=MAX(ixs,0);
            xalpha=xs-ixs;
            xalpha0=1.0-xalpha;

            for (izo=0;izo<gd->nzo;izo++) {
                  zo=gd->fzo+izo*gd->dzo;
    	          zs=(zo-gd->fzs)/gd->dzs;
                  izs=(int)zs;
	          izs=MIN(izs,gd->nzs-2);
	          izs=MAX(izs,0);
                  zgamma=zs-izs;
                  zgamma0=1.0-zgamma;

                  for (iys=0;iys<gd->nys;iys++) {
		        for (iyt=0;iyt<gd->nyt;iyt++) {
			      for (ixt=0;ixt<gd->nxt;ixt++) {
                                    if (ttab[0][iys][ixs  ][izs][iyt][ixt]>=InfTime ||
					ttab[0][iys][ixs+1][izs][iyt][ixt]>=InfTime ||
					ttab[0][iys][ixs  ][izs+1][iyt][ixt]>=InfTime ||
					ttab[0][iys][ixs+1][izs+1][iyt][ixt]>=InfTime) {
				          tbuf[0][iys][ixo][izo][iyt][ixt]=InfByte2;
                                          cbuf[iys][ixo][izo][iyt][ixt]=0;
					  rbuf[iys][ixo][izo][iyt][ixt]=InfByte1;
				    } else {
                                          t=zgamma0*(
                                              xalpha0*ttab[0][iys][ixs  ][izs][iyt][ixt]+
                                              xalpha *ttab[0][iys][ixs+1][izs][iyt][ixt])+
                                            zgamma *(
					      xalpha0*ttab[0][iys][ixs  ][izs+1][iyt][ixt]+
                                              xalpha *ttab[0][iys][ixs+1][izs+1][iyt][ixt]);
                                          tbuf[0][iys][ixo][izo][iyt][ixt]=
					        (unsigned short)((t/InfTime)*InfByte2);

					  if (gd->crfp==NULL) continue;
					  c=zgamma0*(
                                              xalpha0*ctab[iys][ixs  ][izs][iyt][ixt]+
                                              xalpha *ctab[iys][ixs+1][izs][iyt][ixt])+
                                            zgamma *(
					      xalpha0*ctab[iys][ixs  ][izs+1][iyt][ixt]+
                                              xalpha *ctab[iys][ixs+1][izs+1][iyt][ixt]);
					  c=MIN(c,1.0);
					  cbuf[iys][ixo][izo][iyt][ixt]=
					        (unsigned char)(c*InfByte1);

					  r=zgamma0*(
                                              xalpha0*rtab[iys][ixs  ][izs][iyt][ixt]+
                                              xalpha *rtab[iys][ixs+1][izs][iyt][ixt])+
                                            zgamma *(
					      xalpha0*rtab[iys][ixs  ][izs+1][iyt][ixt]+
                                              xalpha *rtab[iys][ixs+1][izs+1][iyt][ixt]);
					  r=MIN(r,InfDistance);
					  rbuf[iys][ixo][izo][iyt][ixt]=
					      (unsigned char)((r/InfDistance)*InfByte1);
				    }
				    for (im=1;im<gd->multit;im++) {
                                        if (ttab[im][iys][ixs  ][izs  ][iyt][ixt]>=InfTime ||
				            ttab[im][iys][ixs+1][izs  ][iyt][ixt]>=InfTime ||
				            ttab[im][iys][ixs  ][izs+1][iyt][ixt]>=InfTime ||
				            ttab[im][iys][ixs+1][izs+1][iyt][ixt]>=InfTime) 
				              tbuf[im][iys][ixo][izo][iyt][ixt]=InfByte2;
				        else {
                                            t=zgamma0*(
                                            xalpha0*ttab[im][iys][ixs  ][izs][iyt][ixt]+
                                            xalpha *ttab[im][iys][ixs+1][izs][iyt][ixt])+
                                              zgamma *(
			                    xalpha0*ttab[im][iys][ixs  ][izs+1][iyt][ixt]+
                                            xalpha *ttab[im][iys][ixs+1][izs+1][iyt][ixt]);
                                            tbuf[im][iys][ixo][izo][iyt][ixt]=
			      	            (unsigned short)((t/InfTime)*InfByte2);
					  }
				    }
			      }
                        }
                  }
	    }
      }
}

void filt( float *trace, struct GD *gd, float *trf)
/**********************************************************************
filt - Low-pass filter, integration and phase shift for input data
***********************************************************************
Function prototype:
void filt( float *trace, struct GD *gd, float *trf);
***********************************************************************
Inputs:
float *trace single seismic trace
struct *gd   grid information

outputs:
float *trf  filtered and phase-shifted seismic trace
***********************************************************************
Author: CWP: Zhaobo Meng, Sept 1997
***********************************************************************/
{
      static int nfft=0; 
      static int itaper, nw, nwf;
      static float *taper, *amp, *ampi, dw;

      int  it, iw, itemp;
      float temp, ftaper, const2, *rt;
      complex *ct;

      gd->fmax *= 2.0*PI;
      ftaper = 0.1*gd->fmax;
      const2 = sqrt(2.0);

      /***************************************************************
      this will only be done once, nice! All the static values are set
      ****************************************************************/
      if (nfft==0) {
            /* Set up FFT parameters */
            nfft = npfaro(gd->nt+gd->ntrimax, 2 * (gd->nt+gd->ntrimax));
            if (nfft >= SU_NFLTS || nfft >= 720720)
                  err("Padded nt=%d -- too big", nfft);

            nw = nfft/2 + 1;
            dw = 2.0*PI/(nfft*gd->dt);

            itaper = 0.5+ftaper/dw;

            taper = ealloc1float(2*itaper+1);

            for (iw=-itaper; iw<=itaper; ++iw){
                  temp = (float)iw/(1.0+itaper);
                  taper[iw+itaper] = (1-temp)*(1-temp)*(temp+2)/4.0;
            }

            nwf = 0.5+gd->fmax/dw;
            if (nwf>nw-itaper-1) nwf = nw-itaper-1;
            amp = ealloc1float(nwf+itaper+1);
            ampi = ealloc1float(nwf+itaper+1);
            amp[0] = ampi[0] = 0.;
            for (iw=1; iw<=nwf+itaper; ++iw){
                  amp[iw] = sqrt(dw*iw)/nfft;
                  ampi[iw] = 0.5/(1-cos(iw*dw*gd->dt));
            }
      }

      /* Allocate fft arrays */
      rt   = ealloc1float(nfft);
      ct   = ealloc1complex(nw);

      memcpy(rt, trace, gd->nt*FSIZE);
      memset((void *) (rt + gd->nt), 0, (nfft- gd->nt)*FSIZE);
      pfarc(1, nfft, rt, ct);

      for(iw=nwf-itaper;iw<=nwf+itaper;++iw){
            itemp = iw-(nwf-itaper);
            ct[iw].r = taper[itemp]*ct[iw].r;
            ct[iw].i = taper[itemp]*ct[iw].i;
      }
      for(iw=nwf+itaper+1;iw<nw;++iw){
            ct[iw].r = 0.;
            ct[iw].i = 0.;
      }

      if(!gd->ls){
            for(iw=0; iw<=nwf+itaper; ++iw){
                  /* phase shifts PI/4    */
                  temp = (ct[iw].r-ct[iw].i)*amp[iw]*const2;
                  ct[iw].i = (ct[iw].r+ct[iw].i)*amp[iw]*const2;
                  ct[iw].r = temp;
            }
      } else {
            for(iw=0; iw<=nwf+itaper; ++iw){
                  ct[iw].i = ct[iw].i*amp[iw];
                  ct[iw].r = ct[iw].r*amp[iw];
            }
      }
      pfacr(-1, nfft, ct, rt);

      /* Load traces back in */
      for (it=0; it<gd->nt; ++it) trace[it] = rt[it];

      /* Integrate traces   */
      for(iw=0; iw<=nwf+itaper; ++iw){
             ct[iw].i = ct[iw].i*ampi[iw];
             ct[iw].r = ct[iw].r*ampi[iw];
      }
      pfacr(-1, nfft, ct, rt);
      for (it=0; it<gd->ntrimax; ++it)  trf[it] = rt[nfft-gd->ntrimax+it];
      for (it=0; it<gd->nt + gd->ntrimax; ++it)  trf[it+gd->ntrimax] = rt[it];

      free1float(rt);
      free1complex(ct);
}
